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We investigate the interplay between ferromagnetic and superconducting order at the edge of a quantum 
spin Hall insulator. Using complementary analytical and self-consistent numerical approaches, we study a 
ferromagnetic Josephson junction and show how the direct coupling between magnetism and the supercon- 
ducting (7(1) phase gives rise to several unusual phenomena which distinguishes the present system from its 
non-topological equivalent. In particular, we demonstrate how the anomalous current-phase relation triggers 
supercurrent-induced magnetization dynamics and also study the spatial localization of the Majorana fermions 
appearing in the junction. 
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The class of materials known as topological insulators 
(TIs), or quantum spin Hall insulators (QSHIs) in two di- 
mensions, represents a new quantum-state of matter [T). Its 
hallmark is the appearance of topologically protected spin- 
polarized edge-states which are robust against time-reversal 
symmetry invariant perturbations. Recently, it has been re- 
alized [2, 3 1 that such systems can host Majorana fermions 
which obey non-Abelian statistics. Thus, the study of TIs is 
presently generating much interest, both from a fundamen- 
tal physical point of view [4] and in terms of future applica- 
tions in fault-tolerant topological quantum computation J3JI2]. 
An intriguing question relates to how different types of long- 
range orders may accommodate to the unusual electronic en- 
vironment in a TL In particular, several works have recently 
focused on how proximity-induced superconductivity is man- 
ifested in TIs J2] [3] [H HTOl . By adding ferromagnetic corre- 
lations via a proximate material with the desired magnetism, 
it has been realized that it is possible to generate Majorana 
fermions in such TI hybrid structures J2] [3). Moreover, it 
was predicted that an anomalous current-phase relation (CPR) 
would appear for a supercurrent flowing through a topologi- 
cal Josephson junction ifTTIl . Interestingly, it follows from this 
fact that a direct coupling between the magnetic and super- 
conducting order exists which can be manipulated simply by 
current-biasing the system. 

In this Rapid Communication, we address a series of un- 
usual phenomena that arise at the edge of a QSHI due to the 
interplay between ferromagnetism and superconductivity that 
have no counterpart in a non-topological equivalent of such a 
system. We employ both an analytical continuum and a self- 
consistent microscopic numerical model, the latter account- 
ing for non-ideal effects such as the partial depletion of the 
superconducting order parameter near the interface regions. 
First, we demonstrate that the predicted anomalous CPR in 
a superconductor | ferromagnet | superconductor (S|F|S) QSHI 
Josephson junction, see Fig.[TJa), is robust towards non-ideal 
effects in the system, but that it still features different behavior 
compared to the prediction of analytical models. Second, we 
demonstrate that this anomalous CPR can be experimentally 



observed and manipulated by means of supercurrent-induced 
magnetization dynamics. This is done by numerically solving 
the full Landau-Lifshitz-Gilbert equation to obtain the time- 
dependence of the magnetization under the influence of an 
AC Josephson effect. Finally, we study the localization of 
Majorana fermions appearing at the S|F interface regions, a 
property crucial with regards to actual experimental schemes 
intending to exploit their non-Abelian properties. 




FIG. 1: (Color online) (a): Schematic of proposed experimental 
setup. Two superconducting contacts (SC) induces superconductivity 
in the S regions of the QSHI whereas the ferromagnetic insulator 
(FM) induces an effective magnetic field in the F region, (b): Kane- 
Mele model for the QSHI with a zigzag edge. For details, see text. 

For the self-consistent microscopic treatment of a QSHI 
S|F|S junction we consider a minimal microscopic model 
H = Hq + Hs + Hp defined on the honeycomb lat- 
tice. Here H = -tJ2(i,j), a a L a j a + Mi Ei, a a \a a i<* + 

iX J2((i j)) l 'i] a \u CT af3 a jl3 i s tne Kane-Mele model for a QSHI 
lfl2ll with a the fermion operator, a, (3 the spin indices, or the 
Pauli matrices, and denoting nearest neighbors 

and next nearest neighbors respectively, and i/y = +1 (— 1) 
if the electron makes a left (right) turn to get to the second 
bond, t = 1 is the nearest neighbor hopping amplitude and 
we also fix the spin-orbit coupling A = 0.3 which gives a 
bulk band gap of 1. Using a unit cell length of c = 1 we get 
a linear dispersion coefficient hvp = 1.15 for the Dirac-like 
edge state which we choose to be located along the zigzag 
edge, see Fig. [TJb) for details. We set the chemical poten- 
tial Hi = /is = 0.3 in the S regions, assumed large due to 
induced doping from the external superconducting contacts, 
whereas we have studied F regions with (j,p = 0,0.1. The 
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superconductivity in the S regions of the QSHI is modeled 
using an attractive Hubbard-?/ pairing potential which is in- 
duced by the external superconducting contacts [ 10 1 such that 
Hs = J2i ^i a !-f a tT a 4 a 4> wnere Ui = U in the S regions but 
zero otherwise. Treating Hs in mean-field theory, the self- 
consistency condition for the superconducting order parame- 
ter A reads Aj = — f7j(a^a^). In general, we choose U such 
that the superconducting coherence length £ = Hvf/Aq in 
the interior of the S regions along the edge satisfies £ ~ L, 
where L is the junction length. Finally, Hp models the in- 
fluence of the external ferromagnet insulator on the QSHI by 
Hp = J2i a ta h i ' Vapaip, where h t = (h x ,h y ,h z ) is the 
(constant) induced magnetic field in the F region of the QSHI, 
ignoring any orbital contribution. Further, we let both the S 
and F regions extend well into the bulk of the QSHI and we 
assume sharp S|F interfaces. We solve the above model self- 
consistently by first diagonalizing the Hamiltonian with the 
order parameter Aj = Ao in the S regions. Aj is then re- 
calculated using the self-consistency condition. By repeating 
this process until Aj does not change between two subsequent 
iterations, we achieve self-consistency for the superconduct- 
ing state and thus capture the inverse proximity effect (IPE), 
i.e. the decrease of A, or equivalently loss of Cooper pairs, on 
the S side of the junction. By applying a difference Acf) in the 
U(l) phase between the order parameters in the two S regions 
we can also calculate the Josephson supercurrent I through 
the junction using the continuity equation for the charge cur- 
rent ifTJl . However, the actual phase drop (f> across the junc- 
tion itself, i.e. across F, necessarily satisfies <j> < A<fi since a 
finite current will always cause a finite phase drop even in the 
S regions. We ensure full self-consistency by only fixing the 
phase in the outermost regions of the S regions. For further 
details on the self-consistency method see Refs. ifTJlfT?! . 

Our non-self-consistent treatment of a QSHI S|F|S junction 
relies on fixing A = Ao in the S regions and using the contin- 
uum field-theoretic model for the QSHI edge state introduced 
by Fu and Kane [3|. The corresponding Hamiltonian is writ- 
ten as: 



H = 



cr ■ (vpk + h) 
-KA)* 



ia y A 



a* ■ (vpk — h) 



(1) 



According to the geometry in Fig. [T] k = k z z. For an arbi- 
trary magnetization direction in the F region, an analytical ex- 
pression for the Andreev-bound states in a short, L < £, S|F|S 
junction can be obtained by constructing the scattering eigen- 
states and matching them at each S|F interface. Performing 
this calculation along the QSHI edge, we obtain the bound- 
state energies e = ±eo, where Bq = AqaAD cos(0/2 — h z L) 
and D = [1 + sinh 2 ^)]" 1 and k = {h 2 x + hi) 1 / 2 . When 
h y = h z = 0, this reduces to the result of Ref. J3J. Interest- 
ingly, the z-component of the exchange field (parallel with the 
junction) effectively renormalizes the superconducting phase- 
difference 1111 . This occurs due to the linear Dirac energy- 
momentum dispersion of the normal-state Hamiltonian for the 
topological edge-state, which effectively makes the exchange 
field enter as a vector potential. The periodic dependence 



of the current on the exchange field is reminiscent from the 
physics of S|F|S systems in conventional metals, with the im- 
portant distinction that these do not depend on the magnetiza- 
tion orientation 1211 . We shall later demonstrate that this fact 
can be both experimentally observed and exploited for practi- 
cal purposes via magnetization dynamics induced by a current 
bias along the edge of the topological insulator. Prior to this 
we will report on the consequences of self-consistency for the 
CPR. 

CPR. The CPR is one of the most important properties 
of a Josephson junction. The non-self-consistent continuum 
model result may be obtained from the above expression for 
the Andreev-bound state e, and gives that I oc Ao sin(<^>/2 — 
h z L)sgn[cos((f>/2 — h z L)] as T — > 0, i.e. the maximum cur- 
rent is reached at the critical phase <p c — > ir for fields along 
the x or y-direction, whereas fields along the z-direction shifts 
this value. This shift is precisely due to the fact that the ex- 
change field enters as a vector potential and thus gives rise to 
an effective flux along the edge. We note in passing that in the 
non-thermalized regime, it is possible to obtain a Air -periodic 
CPR on a topological surface (3). 

The self-consistent microscopic results for non-zero h x and 
h z is shown Figs.^a) and (b), respectively. A finite /i^-field 
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FIG. 2: (Color online) The supercurrent I/Io with Jo = eAo/fi 
and { = 13 as function of phase drop (j> across the junction for (a): 
h x = 0.3, fi F = 0.1, and L = 12 (black) and L = 24 (red) with 
dashed lines showing results after one self-consistency step only, and 
(b): h x = (black,x), 0.01 (red,o), -0.01 (green,D), fi F = 0, and 
L — 12. Inset shows the IPE at the S|F interface by displaying 
A/A for L = 12, £ = 13 (black), £ = 5 (dashed) and L = 24, 
£ = 13 (red). 

causes Majorana modes to appear at the S|F interfaces J3J and 
we show in Fig. [2ja) results for both L — 0.9£ and L = 1.8£. 
We clearly see that C is shifted down toward 7r/2 for both 
cases. We attribute this shift to a significant amount of IPE, as 
seen in the inset. In fact, processes influencing the S regions 
are known to shift cf) c below the standard value of 7r/2 in con- 
ventional S | N | S junctions [ 1 5 1 . Here we see a similar shift but 
starting from the non-self-consistent solution at (f> c — > tt. To 
support this conclusion we show in the dashed lines the CPR 
after the first step of the self-consistency loop, starting from 
the analytical A = Ao step-function solution. Already af- 
ter one iteration we see that <j) c is significantly reduced while 
at the same time a sizable IPE has appeared. We can quan- 
tify the IPE by fitting A close to the interface to the form 
A = Aq[1 — exp(— z/kipe)]. As expected, the decay length 
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Kipe is linearly proportional to £, although the linear coeffi- 
cient is somewhat dependent on both L and h x , but not on [ip. 
Thus, by varying £ we can significantly modify the IPE and for 
very short £ we indeed see C increasing, but never fully to tt. 
On the other hand, we only achieve C < 7r/2 due to current 
depairing fi6l . which is only present in artificially short junc- 
tions. Therefore, for any experimental junction, where £ is 
significantly larger than we can model, we predict C « 7r/2. 

Figure |2|b) shows the results for finite h z -fields. There are 
no Majorana fermions appearing for this field direction, but 
still, a h z -field has interesting consequences for the CPR. As 
predicted analytically, we see a linear shift in C with h z and 
L, however, the linear coefficient is ~ 1.5. Moreover, we 
again see C w tt/2, and not tt, in the limit of zero field. The 
self-consistent CPR with a /i z -field can thus, up to numerical 
accuracy, be written as / = I c sin(0 — 1.5h z L). We also see 
a shift in I c with h z (but not L), such that I c = I c (h z — 
0) + 0.lh z . Such a dependence is not seen in the analytical 
solution but can be traced back to a changed IPE. In fact, when 
the IPE increases, we see a concurrent drop of both I c and 
4> c , which can also explain the larger linear coefficient for the 
field-induced phase shift than found analytically. 

Magnetization dynamics. The direct coupling between the 
superconducting phase difference and the exchange-field of- 
fers an interesting opportunity: supercurrent-induced mag- 
netization dynamics. By biasing the system with a super- 
current, one would expect to see a time-evolution of the 
magnetic order parameter. To address this issue, we make 
use of the Landau-Lifshitz-Gilbert (LLG) equation ifTTl (the 
modification of the LLG equation on a topological surface 
in the presence of a single-particle charge-current was re- 
cently discussed in Ref. [18]): d t M = -7M x H eS + 
(a/Mo)(M x dtM), where M is the magnetization vec- 
tor, Mo = 7 is the gyromagnetic ratio, and a is the 
Gilbert damping constant. The components of the effective 
field are obtained from the free energy T of the system via: 
(-ffeff)j = — A , where V is the volume of the F region. 
In what follows, we investigate how an equilibrium supercur- 
rent can induce strong magnetization dynamics in the present 
system. Thus, we focus on the effect of the Andreev-bound 
states, the rationale being that it is precisely this contribu- 
tion that constitutes the Josephson-current induced magneti- 
zation dynamics, and assume that other contributions such 
as anisotropy fields can be neglected. We underline the fact 
that the analytical expression for the bound-states captures 
the field-dependence exhibited by the effective superconduct- 
ing phase difference, as also confirmed by our self-consistent 
results, which is the key to the supercurrent-induced mag- 
netization dynamics. Obtaining the free energy from T = 
— ^ln[2 cosh(/3eo/2)], one may rewrite the LLG equation in 
the convenient dimensionless form: d T n = —n x H e a + 
a(n x d T n), where we have defined n = (n x ,n y ,n z ) — 
M/\M\, t = ujjt where ojj is the Josephson fre- 
quency, and finally: ("H e ff)j = C\Uj tanh(/3e /2) cos(0/2 — 
hLn z ), j = x,y, (7i eS ) z = C 2 tanh(/3e /2) sin (0/2 - 
hLn z ). Here, d = -A h 2 Lsmh(kL)/(2kD 5 / 2 ujT) and 




FIG. 3: (Color online) Time-evolution of the magnetization compo- 
nents for (a) a — 0.01 and (b) a — 0.1. Thick dashed line: n z , thick 
full line: n y , thin full line: n x . The top panels show the magnetiza- 
tion dynamics on a larger time-scale. 



C 2 = A v r DhL/(2cjjr) where T = M V/(2 7 ). The in- 
verse temperature is (3, and we set T/T c = 0.2. 

We have solved the above LLG equation numerically to 
simulate the magnetization dynamics induced by a Josephson 
current. It should be noted that a Josephson current does not 
induce magnetization dynamics in a conventional S|F|S junc- 
tion, whereas this feature appears in the present context due 
to the Dirac-like dispersion with momentum-spin locking for 
the topological edge states. Let us estimate the magnitude of 
the parameters entering in the expressions for C\ and C2. The 
proximity induced superconducting gap is taken to be small, 
A = 0.1 meV, and we also make a conservative assumption 
for the exchange field, h = 5 meV (Ref. 11221 ). The length L 
of the junction can be varied experimentally, but is here fixed 
to L = 400 nm. We use a typical value for the Josephson 
frequency, u>j = 1 GHz. The parameter a is phenomenolog- 
ical, and we will contrast the small damping scenario a -C 1 
with strong damping a ~ 0.1. We also treat V as a free pa- 
rameter to vary, since the exact values of Mq and V will differ 
depending on the specific materials and dimensions chosen. 
When r> 1, the effective field in the junction is reduced and 
the magnetization dynamics vanishes. Therefore, we here set 
r = 0.2 and compare in Fig. |3]the case of weak (a = 0.01) 
and strong (a = 0.1) damping. The initial condition is chosen 
to have the magnetization || y, i.e. in the plane of the topolog- 
ical insulator. In the weak damping regime, the magnetization 
oscillates unabatedly as t increases (only the main periodic 
pattern is shown in Fig. [3] for clarity). For larger damp- 
ing, the dynamics becomes more interesting. In this case, 
it is seen how the z-component becomes dominant and how 
the oscillation of the other components slowly decay as r in- 
creases. This can be understood from the fact that the magne- 
tization will saturate on a shorter time-scale when a <C 1 is no 
longer satisfied. In fact, the right panel of Fig. |3]suggests that 
supercurrent-induced magnetization switching might be pos- 
sible: the magnetization starts out || y, but is seen to converge 
toward || z as r increases. An important point with respect to 
the presence of an exchange field in the z-direction is that the 
ground-state superconducting phase difference 4> no longer 
is 0o = or 0o = if, as in conventional metallic systems. 
Instead, the ground-state phase is now given by 0o = 2h z L. 
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This is similar to the scenario of a noncentrosymmetric fer- 
romagnet in a Josephson junction fl9\ , although the origin of 
the 0o-state here is completely different. We also note that the 
same magnetization dynamics effect would be present on the 
surface of three-dimensional (3D) topological insulators, pos- 
sibly exerting stronger dynamics in magnitude in that case due 
to the larger number of Andreev-bound states carrying cur- 
rent. Finally, it should be noted that we have assumed that the 
magnetic anisotropy energy is small enough to be negligible, 
which places restrictions on which type of materials that may 
be used for this purpose. This fact notwithstanding, the prin- 
ciple of supercurrent-induced magnetization dynamics would 
remain possible due to the renormalized phase-difference. 

Majorana fermions. Our self-consistent solution also of- 
fers access to various properties of the Majorana fermions 
which appear at the S|F interfaces at finite h x -fields. Here 
we concentrate on their spatial extension (i.e. the degree of 
localization to the S|F interface) at <fr = 0. This property 
is important for any experimental scheme aimed at exploit- 
ing their non-Abelian properties through braiding. Note that 
we are here not concerned with their spread into the QSHI as 
that is set by the width of the edge state and is therefore very 
narrow. In a S|F|S junction, the two Majorana fermions ap- 
pear as two in-gap, near-zero energy solutions with eigenen- 
ergies ±e and eigenvectors T±, where particle-hole symme- 
try dictates T + = T~[_ = T . A finite L gives rise to a 
finite mixing between the two Majorana modes which, at 
= 0, leaves a non-zero e J3]. Also, the two eigenvectors 
contain an equal mixture of the two Majorana modes living 
at the two separate interfaces. But, by defining the vectors 
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7i = l/2(e^r + e- l »rl) and l2 = i/2(e l0 T 
with a variable 9, we can always isolate each Majorana mode, 
which then also explicitly displays their 7^ = 7 Majorana na- 
ture. By studying the (normalized) density \ji\ 2 , we find that 
the spatial extent of the Majorana mode on the S side of the 
SF interface is only a function of £, whereas on the F side it is 
only a function of h x . Notably, we find no significant depen- 
dence on jjLp or L. Figure |4|a) shows the decay of the Majo- 




FIG. 4: (Color online) Majorana mode density (normalized to 1) 
along the QSHI edge as function of the distance z (unit cells) from 
the S|F interface for L — 12, fiF = 0, and (a): h x = 0.3 with 
£ = 3, 5, 7, 10 (increasing decay length in S) and (b): £ = 13 with 
h x = 0.3,0.2,0.15,0.08 (increasing decay length in F). Inset (a): 
Crosses show km.s vs. £ for L = 8, 12, 24. Line is km,s = 0.6£. 
Inset (b): km,f vs. h x . Line is only a guide to the eye. 

rana mode into the S region for a few representative values of 



£. Note that the slight difference in the curves on the F side for 
different £ is solely due to the normalization constraint, giv- 
ing different peak heights at the interface. The decay on the S 
side can to good accuracy be approximated as an exponential 
decay with a decay length km.s- We find that km.s — 0.6£ for 
all junctions we have studied (see inset). This decay length 
is similar to ripe and one could speculate that the Majorana 
mode extends into the IPE depression region of the supercon- 
ducting order parameter. However, kipe also has a moderate 
dependence on L and h x 1231 which km,s does not. Thus the 
Majorana decay into the S region is not solely governed by the 
IPE although the two show a strong correlation. In Fig. |4|b) 
we show a similar set of data for the decay of |7^| 2 on the F 
side of the junction for different values of h x . As seen in the 
inset, we find a non-linear dependence on h x for km.f ( or the 
energy gap produced in the edge state due to h x ) although it 
is monotonic. 

In summary, we have considered the interplay between fer- 
romagnetism and superconductivity at the edge of a QSHI. We 
have shown that the supercurrent flowing in a QSHI ferromag- 
netic Josephson junction exhibits several unusual phenomena, 
which makes it distinct from its non-topological equivalent. 
Using both analytical and fully self-consistent numerical cal- 
culations, we have studied an anomalous CPR, supercurrent- 
induced magnetization dynamics, and the appearance of Ma- 
jorana fermions along the QSHI edge. Besides the application 
of Majorana fermions to topological quantum computing, the 
interplay between magnetism and the superconducting 17(1) 
phase in this type of system opens new perspectives related to 
spin-polarized and tunable supercurrents. 

A.M.B.-S. thanks Hans Hansson and Eddy Ardonne for 
valuable discussions. 
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